/*
 * @Description: 
 * @Autor: kaikaima
 * @Date: 2021-04-26 20:38:46
 * @LastEditors: kaikaima
 * @LastEditTime: 2021-05-17 14:14:41
 */
#include <iostream>
#include <cmath>

#ifndef __hhxxttxs_MATRIX__
#define __hhxxttxs_MATRIX__

class Matrix{
    private:
    double* house(double* _x, int _n){
        double* v=new double[_n+1];
        double x_norm=0, sigma=0;
        for(int i=0;i<_n;i++)
            x_norm=((fabs(_x[i])>x_norm)?fabs(_x[i]):x_norm);
        for(int i=0;i<_n;i++){
            _x[i]/=x_norm;
            if(i!=0) sigma+=(_x[i]*_x[i]);
            v[i]=_x[i];
        }
        if(sigma==0){
            v[_n]=0;
        }
        else{
            double alpha=sqrt(_x[0]*_x[0]+sigma);
            if(_x[0]<=0) v[0]=_x[0]-alpha;
            else v[0]=-sigma/(_x[0]+alpha);
            v[_n]=2*v[0]*v[0]/(sigma+v[0]*v[0]);
            for(int i=1;i<_n;i++) v[i]/=v[0];
            v[0]=1;
        }
        return v;
    };
    double* House(double* _b, int _i){
        double* result=new double[M];
        for(int i=0;i<M;i++){
            if(i<_i) result[i]=0;
            else if(i==_i){
                result[i]=_b[i];
                for(int j=_i+1;j<M;j++)
                result[i]+=(A[j*N+_i]*_b[j]);
            }
            else{
                result[i]=_b[_i];
                for(int j=_i+1;j<M;j++)
                result[i]+=(A[j*N+_i]*_b[j]);
                result[i]*=A[i*N+_i];
            }
        }
        return result;
    };
    double* QR_Householder(){
        double* d=new double[N];
        for(int i=0;i<N;i++){
            if(i>=M) break;
            double* _x=new double[M-i];
            for(int j=i;j<M;j++)
                _x[j-i]=A[j*N+i];
            double* v=house(_x,M-i);
            delete[] _x;

            double w[N-i]={0};
            for(int _i=i;_i<N;_i++)
            for(int _j=i;_j<M;_j++)
            w[_i-i]+=(A[_j*N+_i]*v[_j-i]*v[M-i]);
            for(int _i=i;_i<M;_i++)
            for(int _j=i;_j<N;_j++)
            A[_i*N+_j]-=(v[_i-i]*w[_j-i]);

            d[i]=v[M-i];
            for(int j=i+1;j<M;j++)
                A[j*N+i]=v[j-i];
            delete[] v;
        }
        return d;
    };

    protected:
    unsigned N;
    unsigned M;
    double *A;

    public:
    Matrix(int _m, int _n){
        M=_m;
        N=_n;
        A=new double [M*N];
    };
    Matrix(int _n) : Matrix(_n, _n) {};
    ~Matrix(){
        delete[] A;
    };
    void set(double **_a){
        for(int i=0;i<M;i++)
        for(int j=0;j<N;j++)
            A[i*N+j]=_a[i][j];
    };
    double* operator*(double* _a) const{
        double* result=new double [M];
        for(int i=0;i<M;i++){
            result[i]=0;
            for(int j=0;j<N;j++)
                result[i]+=(A[i*N+j]*_a[j]);
        }
        return result;
    };
    void transpose(){
        for (int i=0;i<M*N;i++){
            int next=(i%N)*M+i/N;
            while (next > i) next=(next%N)*M+next/N;
            if (next == i){
                int temp = A[i];
                int pre = (i%M)*N+i/M;
                int cur=i;
                while (pre!=i){
                    A[cur] = A[pre];
                    cur = pre;
                    pre = (cur%M)*N+cur/M;
                }
                A[cur] = temp;
            }
        }
        unsigned tmp=M;
        M=N;
        N=tmp;
    };
    /*double norm() const{
        double result=0;
        double w[N],v[N],z[N];
        double x[N]={1,0};
        while(true){
            for(int i=0;i<N;i++){
                w[i]=0;
                for(int j=0;j<N;j++)
                    w[i]+=(A[i*N+j]*x[j]);
            }
            for(int i=0;i<N;i++){
                if(w[i]==0) v[i]=0;
                else if(w[i]>0) v[i]=1;
                else v[i]=-1;
            }
            for(int i=0;i<N;i++){
                z[i]=0;
                for(int j=0;j<N;j++)
                z[i]+=(A[i*N+j]*v[j]);
            }
            //get ||z||∞ zTx
            int zn=0;
            double zx=0;
            for(int i=1;i<N;i++) if(fabs(z[i])>fabs(z[zn])) zn=i;
            for(int i=0;i<N;i++) zx+=(z[i]*x[i]);

            if(fabs(z[zn])<=zx){
                for(int i=0;i<N;i++)
                result+=fabs(w[i]);
                break;
            }
            else{
                for(int i=0;i<N;i++)
                    x[i]=(i==zn?1:0);
            }
        }
        return result;
    };
    */
    void print(){
        for(int i=0;i<M;i++){
            for(int j=0;j<N;j++)
            std::cout<<A[i*N+j]<<" ";
            std::cout<<std::endl;
        }
        std::cout<<std::endl;
    };
    double* Least_squares(double* _b){
        double* b=QR_Householder();
        for(int i=0;i<N;i++){
            double* v=House(_b,i);
            for(int j=0;j<M;j++)
                _b[j]-=(v[j]*b[i]);
            delete[] v;
        }
        delete[] b;
        for(int i=N-1;i>0;i--){
            _b[i]/=A[i*N+i];
            for(int j=0;j<i;j++)
                _b[j]-=(_b[i]*A[j*N+i]);
        }
        _b[0]/=A[0];
        double* result=new double[N];
        for(int i=0;i<N;i++) result[i]=_b[i];
        delete[] _b;
        return result;
    };
};

#else
//DO NOTHING!
#endif